// STL
#include <cmath>                         // std::sqrt()
#include <limits>                        // std::numeric_limits<>
#include <stdexcept>                     // std::runtime_error()
#include <vector>                        // std::vector<>

// functions++
#include "functionsxx/integration.hpp"   // functionsxx::integration::curve_distances()
#include "functionsxx/interpolation.hpp" // functionsxx::interpolation::spline(), ::splint()

// stats++
#include "statsxx/postprocess/ROC.hpp"   // ROC_pt


//
// DESC: Map (discrete) ROC scores to an ROC curve.
//
// -----
//
// NOTE: This subroutine maps the (discrete) scores to the ROC curve, and then does a spline interpolation along all points.
//
// -----
//
// NOTE: This subroutine, and the associated that are called, assume that the ROC curve can be constructed as piecewise linear.
//
inline std::vector<ROC_pt> map_scores(
                                      const std::vector<ROC_pt> &ROC_scores,
                                      // -----
                                      std::vector<ROC_pt>        ROC
                                      )
{
//    std::vector<ROC_pt> ROC;

    //=========================================================
    // CALCULATE DISTANCES ALONG ROC CURVE
    //=========================================================

    std::vector<double> x;
    std::vector<double> y;

    for( auto i = 0; i < ROC.size(); ++i )
    {
        x.push_back(ROC[i].FPR);
        y.push_back(ROC[i].TPR);
    }

    std::vector<double> d = functionsxx::integration::curve_distances(
                                                                      x,
                                                                      y
                                                                      );

    //=========================================================
    // MAP (DISCRETE) ROC SCORES (ROC_scores) TO CORRESPONDING POINTS ON THE ROC CURVE
    //=========================================================

    std::vector<double> d_scores;
    std::vector<double> scores;

    for( auto i = 0; i < ROC_scores.size(); ++i )
    {
        // FIND CLOSEST POINT ON ROC CURVE TO POINT
        bool found = false;

        int    idx1;
        int    idx2;

        double FPR;
        double TPR;

        double d2min = std::numeric_limits<double>::max();

        for( auto j = 1; j < ROC.size(); ++j )
        {
            // GET THE EQUATION OF THE LINE PASSING THROUGH [i-1] AND [i]
            //
            // NOTE:
            //
            //     y = mx + b
            //
            // NOTE: See: http://mathworld.wolfram.com/Two-PointForm.html

            double m = (ROC[j].TPR - ROC[j-1].TPR)/(ROC[j].FPR - ROC[j-1].FPR);

            double b = ROC[j-1].TPR - m*ROC[j-1].FPR;

            // FIND THE POINT ON THIS LINE CLOSEST TO (SCORE) POINT
            //
            // NOTE:
            //
            //     ax + by + c = 0
            //
            // NOTE: See: https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line

            double _a = m;
            double _b = -1.;
            double _c = b;

            double _FPR = (_b*(_b*ROC_scores[i].FPR - _a*ROC_scores[i].TPR) - _a*_c)/(_a*_a + _b*_b);
            double _TPR = (_a*(-_b*ROC_scores[i].FPR + _a*ROC_scores[i].TPR) - _b*_c)/(_a*_a + _b*_b);

            // CHECK WHETHER THIS POINT IS IN THE ROC INTERVAL
            //
            // NOTE: The FPR always increases from left to right, ...
            // NOTE: ... whereas (for a very improper curve) TPR may decrease.
            //
            // NOTE: The <=/>= on both ends allows this to capture the end points, and should not cause a problem if the point lies on a point of a curve.
            if( (_FPR >= ROC[j-1].FPR) && (_FPR <= ROC[j].FPR) )
            {
                found = true;

                // ... AND IF SO, CALCULATE THE DISTANCE FROM (SCORE) POINT TO LINE
                double d2 = (ROC_scores[i].TPR - _TPR)*(ROC_scores[i].TPR - _TPR) + (ROC_scores[i].FPR - _FPR)*(ROC_scores[i].FPR - _FPR);

                if( d2 < d2min )
                {
                    idx1 = j - 1;
                    idx2 = j;

                    FPR = _FPR;
                    TPR = _TPR;

                    d2min = d2;
                }
            }
        }

        if( !found )
        {
            throw std::runtime_error("error in map_scores(): could not map score to ROC curve");
        }

        // CALCULATE DISTANCE OF THE MINIMUM POINT (FOUND ABOVE) ALONG THE ROC CURVE
        double _d = d[idx1] + std::sqrt((FPR - ROC[idx1].FPR)*(FPR - ROC[idx1].FPR) + (TPR - ROC[idx1].TPR)*(TPR - ROC[idx1].TPR));

        d_scores.push_back(_d);
        scores.push_back(ROC_scores[i].score);
    }

    //=========================================================
    // SPLINE INTERPOLATE (DISCRETE) SCORES TO ROC CURVE
    //=========================================================
    //
    // NOTE: A natural spline interpolation is used.

    std::vector<bool>   natural(2, true);
    std::vector<double> yp(2); // unused for natural spline
    std::vector<double> y2 = functionsxx::interpolation::spline(
                                                                d_scores,
                                                                scores,
                                                                natural,
                                                                yp
                                                                );

    for( auto i = 0; i < ROC.size(); ++i )
    {
        ROC[i].score = functionsxx::interpolation::splint(
                                                          d_scores,
                                                          scores,
                                                          y2,
                                                          d[i]
                                                          );
    }

    //=========================================================
    // RETURN
    //=========================================================

    return ROC;
}
